
function G = grad_cfs_moments(theta,X_PVEM,X_PRI,ES0,ES1,ES2,DES0,DES1,DES2,DPW0_PVEM,DPW0_PRI)

P1 = exp(X_PVEM * theta + ES1 - ES0);
P2 = exp(X_PRI * theta + ES2 - ES0);
P0 = 1 + P1 + P2;
P1 = P1 ./ P0;
P2 = P2 ./ P0;

DP1 = theta(end) * DPW0_PVEM + DES1 - DES0;
DP2 = theta(end) * DPW0_PRI + DES2 - DES0;

H = HessLL_cfs(theta,0,X_PVEM,X_PRI,ES0,ES1,ES2);

G = - X_PVEM' * diag(P1.*(1-P1)) * DP1 + ...
    X_PVEM' * diag(P1.*P2) * DP2 + ...
    X_PRI' * diag(P1.*P2) * DP1 + ...
    - X_PRI' * diag(P2.*(1-P2)) * DP2;

G=[H,-G];



